/*
 * System.cpp
 *
 * Created on: 23 May 2011
 *   Author: allan
 */

#include "System.h"

// C++ includes
#include <cmath>
#include <iostream>
#include <set>
using namespace std;

System::System() :
aqueousPhase(0), gaseousPhase(0)
{}

System::System(const System& copy)
{
	
}

System::~System()
{
	for(unsigned iPhase = 0; iPhase < GetNumOfPhases(); ++iPhase)
		delete phases[iPhase];	
}

void System::AddPhase(const AqueousPhase& aqueousPhase)
{
	this->aqueousPhase = new AqueousPhase(aqueousPhase);
	
	phases.push_back(this->aqueousPhase);
	
	activitiesOfPhases.push_back(vector<double>(aqueousPhase.GetNumOfSpecies()));
	
	UpdateSpeciesWithNewPhase();
}

void System::AddPhase(const GaseousPhase& gaseousPhase)
{
	this->gaseousPhase = new GaseousPhase(gaseousPhase);
	
	phases.push_back(this->gaseousPhase);
	
	activitiesOfPhases.push_back(vector<double>(gaseousPhase.GetNumOfSpecies()));
	
	UpdateSpeciesWithNewPhase();
}

void System::AddPhase(const MineralPhase& mineralPhase)
{
	mineralPhases.push_back(new MineralPhase(mineralPhase));
	
	phases.push_back(mineralPhases.back());
	
	activitiesOfPhases.push_back(vector<double>(mineralPhase.GetNumOfSpecies()));
	
	UpdateSpeciesWithNewPhase();
}

void System::AddEquilibriumReaction(const EquilibriumReaction& equilibriumReaction)
{
	equilibriumReactions.push_back(equilibriumReaction);
	
	unsigned indexOfReactant = unsigned(-1);
	
	for(unsigned iEquilibrium = 0; iEquilibrium < GetNumOfEquilibriumSpecies(); ++iEquilibrium)
		if(species[idxEquilibriumSpecies[iEquilibrium]].GetName() == equilibriumReaction.GetNameOfReactant())
			indexOfReactant = iEquilibrium;
	
	vector<unsigned> indexesOfProducts;
	
	for(unsigned iProduct = 0; iProduct < equilibriumReaction.GetNumOfProducts(); ++iProduct)
		for(unsigned iPrimary = 0; iPrimary < GetNumOfPrimarySpecies(); ++iPrimary)
			if(species[idxPrimarySpecies[iPrimary]].GetName() == equilibriumReaction.GetNameOfProduct(iProduct))
				indexesOfProducts.push_back(iPrimary);
	
	equilibriumReactions.back().SetIndexOfReactant(indexOfReactant);
	equilibriumReactions.back().SetIndexOfProducts(indexesOfProducts);
}

void System::AddKineticReaction(const MineralReaction& mineralReaction)
{
	kineticReactions.push_back(new MineralReaction(mineralReaction));
	
	unsigned indexOfReactant = unsigned(-1);
	
	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
		if(species[idxKineticSpecies[iKinetic]].GetName() == mineralReaction.GetNameOfReactant())
			indexOfReactant = iKinetic;
	
	vector<unsigned> indexesOfProducts;
	
	for(unsigned iProduct = 0; iProduct < mineralReaction.GetNumOfProducts(); ++iProduct)
		for(unsigned iPrimary = 0; iPrimary < GetNumOfPrimarySpecies(); ++iPrimary)
			if(species[idxPrimarySpecies[iPrimary]].GetName() == mineralReaction.GetNameOfProduct(iProduct))
				indexesOfProducts.push_back(iPrimary);
	
	kineticReactions.back()->SetIndexOfReactant(indexOfReactant);
	kineticReactions.back()->SetIndexOfProducts(indexesOfProducts);
}

void System::SetPrimarySpecies(const vector<string>& namesOfPrimarySpecies)
{
	for(unsigned iter = 0; iter < namesOfPrimarySpecies.size(); ++iter)
		idxPrimarySpecies.push_back(GetIndexOfSpecies(namesOfPrimarySpecies[iter]));
}

void System::SetEquilibriumSpecies(const vector<string>& namesOfEquilibriumSpecies)
{
	for(unsigned iter = 0; iter < namesOfEquilibriumSpecies.size(); ++iter)
		idxEquilibriumSpecies.push_back(GetIndexOfSpecies(namesOfEquilibriumSpecies[iter]));
	
	Qe.resize(idxEquilibriumSpecies.size());
}

void System::SetKineticSpecies(const vector<string>& namesOfKineticSpecies)
{
	for(unsigned iter = 0; iter < namesOfKineticSpecies.size(); ++iter)
		idxKineticSpecies.push_back(GetIndexOfSpecies(namesOfKineticSpecies[iter]));
	
	Qk.resize(idxKineticSpecies.size());
}

void System::SetMolesOfSpecies(unsigned iSpecies, double value)
{
	phases[species[iSpecies].GetIdxPhase()]->SetMolesOfSpecies(species[iSpecies].GetIdxLocal(), value);
}

void System::SetMolesOfSpecies(const vector<double>& m)
{
	for(unsigned iSpecies = 0; iSpecies < GetNumOfSpecies(); ++iSpecies)
		SetMolesOfSpecies(iSpecies, m[iSpecies]);
}

void System::SetMolesOfPrimarySpecies(const vector<double>& mj)
{
	for(unsigned iPrimary = 0; iPrimary < GetNumOfPrimarySpecies(); ++iPrimary)
		SetMolesOfSpecies(idxPrimarySpecies[iPrimary], mj[iPrimary]);
}

void System::SetMolesOfEquilibriumSpecies(const vector<double>& me)
{
	for(unsigned iEquilibrium = 0; iEquilibrium < GetNumOfEquilibriumSpecies(); ++iEquilibrium)
		SetMolesOfSpecies(idxEquilibriumSpecies[iEquilibrium], me[iEquilibrium]);
}

void System::SetMolesOfKineticSpecies(const vector<double>& mk)
{
	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
		SetMolesOfSpecies(idxKineticSpecies[iKinetic], mk[iKinetic]);
}

void System::SetMolesOfPhase(unsigned iPhase, const vector<double>& m)
{
	phases[iPhase]->SetMolesOfSpecies(m);
}

void System::SetTemperature(double T)
{
	this->T = T;
}

void System::SetPressure(double P)
{
	this->P = P;
}

const unsigned System::GetNumOfSpecies() const
{
	return species.size();
}

const unsigned System::GetNumOfPrimarySpecies() const
{
	return idxPrimarySpecies.size();
}

const unsigned System::GetNumOfEquilibriumSpecies() const
{
	return idxEquilibriumSpecies.size();
}

const unsigned System::GetNumOfKineticSpecies() const
{
	return idxKineticSpecies.size();
}

const unsigned System::GetNumOfPhases() const
{
	return phases.size();
}

const unsigned System::GetIndexOfSpecies(string nameOfSpecies) const
{
	for(unsigned iSpecies = 0; iSpecies < species.size(); ++iSpecies)
		if(species[iSpecies].GetName() == nameOfSpecies)
			return iSpecies;
	
	cerr << "Error: There is no species with name " << nameOfSpecies << "." << endl;
	exit(0);
}

const unsigned System::GetIndexOfPhase(string nameOfPhase) const
{
	for(unsigned iPhase = 0; iPhase < phases.size(); ++iPhase)
		if(phases[iPhase]->GetName() == nameOfPhase)
			return iPhase;
	
	cerr << "Error: There is no phase with name " << nameOfPhase << "." << endl;
	exit(0);
}

const unsigned System::GetIndexOfPrimarySpecies(unsigned iPrimary) const
{
	return idxPrimarySpecies[iPrimary];
}

const unsigned System::GetIndexOfEquilibriumSpecies(unsigned iEquilibrium) const
{
	return idxEquilibriumSpecies[iEquilibrium];
}

const unsigned System::GetIndexOfKineticSpecies(unsigned iKinetic) const
{
	return idxKineticSpecies[iKinetic];
}

const Species& System::GetSpecies(unsigned iSpecies) const
{
	return species[iSpecies];
}

const Phase* System::GetPhase(unsigned iPhase) const
{
	return phases[iPhase];
}

const AqueousPhase* System::GetAqueousPhase() const
{
	return aqueousPhase;
}

const GaseousPhase* System::GetGaseousPhase() const
{
	return gaseousPhase;
}

const MineralPhase* System::GetMineralPhase(unsigned iMineralPhase) const
{
	return mineralPhases[iMineralPhase];
}

const EquilibriumReaction& System::GetEquilibriumReaction(unsigned iEquilibrium) const
{
	return equilibriumReactions[iEquilibrium];
}

const KineticReaction& System::GetKineticReaction(unsigned iKinetic) const
{
	return *kineticReactions[iKinetic];
}

const double System::GetMolesOfSpecies(unsigned iSpecies) const
{
	return phases[species[iSpecies].GetIdxPhase()]->GetMolesOfSpecies(species[iSpecies].GetIdxLocal());
}

const double System::GetMolarFractionOfSpecies(unsigned iSpecies) const
{
	return phases[species[iSpecies].GetIdxPhase()]->GetMolarFractionOfSpecies(species[iSpecies].GetIdxLocal());
}

const double System::GetActivityOfSpecies(unsigned iSpecies) const
{
	return activitiesOfPhases[species[iSpecies].GetIdxPhase()][species[iSpecies].GetIdxLocal()];
}

const double System::GetEquilibriumReactionQuotient(unsigned iEquilibrium) const
{
	return Qe[iEquilibrium];
}

const double System::GetKineticReactionQuotient(unsigned iKinetic) const
{
	return Qk[iKinetic];
}

const double System::GetTemperature() const
{
	return T;
}

const double System::GetPressure() const
{
	return P;
}

const vector<double> System::GenerationRatesOfKineticSpecies() const
{
	vector<double> rates(GetNumOfKineticSpecies());
	
	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
	{
		rates[iKinetic] = kineticReactions[iKinetic]->ReactionRate(T, P, *this);
	}
	
	return rates;
}

const vector<double> System::GenerationRatesOfPrimarySpecies(const vector<double>& ratesOfKineticSpecies) const
{
	vector<double> rates(GetNumOfPrimarySpecies(), 0.0);
	
	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
	{
		const KineticReaction& kineticReaction = GetKineticReaction(iKinetic);
		
		for(unsigned iProduct = 0; iProduct < kineticReaction.GetNumOfProducts(); ++iProduct)
		{
			const unsigned iPrimary = kineticReaction.GetIndexOfProduct(iProduct);
			
			const double vkj = kineticReaction.GetStoichiometryOfProduct(iProduct);
			
			rates[iPrimary] -= vkj * ratesOfKineticSpecies[iKinetic];
		}
	}
	
	return rates;
}

void System::UpdateSystem()
{
	UpdateActivities();
	
	UpdateReactionQuotients();
}

void System::UpdateSpeciesWithNewPhase()
{
	// a pointer  to the new phase
	const Phase* newPhase = phases.back();
	
	// the number of species in the phase
	const unsigned numNewSpecies = newPhase->GetNumOfSpecies();
	
	// the index of the new phase in the reactive system
	const unsigned indexOfPhase = phases.size() - 1;
	
	// iterate over all species in the new phase
	for(unsigned iLocalSpecies = 0; iLocalSpecies < numNewSpecies; ++iLocalSpecies)
	{
		Species newSpecies;
		
		newSpecies.SetName(newPhase->GetNameOfSpecies(iLocalSpecies));
		
		newSpecies.SetIdxGlobal(GetNumOfSpecies());
		
		newSpecies.SetIdxLocal(iLocalSpecies);
		
		newSpecies.SetIdxPhase(indexOfPhase);
		
		species.push_back(newSpecies);
	}
}

void System::UpdateActivities()
{
	for(unsigned iPhase = 0; iPhase < GetNumOfPhases(); ++iPhase)
		activitiesOfPhases[iPhase] = phases[iPhase]->Activities(T, P);
}

void System::UpdateReactionQuotients()
{
	for(unsigned iEquilibrium = 0; iEquilibrium < GetNumOfEquilibriumSpecies(); ++iEquilibrium)
		Qe[iEquilibrium] = EquilibriumReactionQuotient(iEquilibrium);
	
	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
		Qk[iKinetic] = KineticReactionQuotient(iKinetic);
}

const double System::EquilibriumReactionQuotient(unsigned iEquilibrium) const // TODO: Modify the design of the class such that there is only one method to compute reaction quotient.
{
	const EquilibriumReaction& reaction = equilibriumReactions[iEquilibrium];
	
	const unsigned r = reaction.GetIndexOfReactant();
	
	const double ar = GetActivityOfSpecies(idxEquilibriumSpecies[r]);
	
	double Qr = 1.0/ar;
	
	for(unsigned iProduct = 0; iProduct < reaction.GetNumOfProducts(); ++iProduct)
	{
		const unsigned j = reaction.GetIndexOfProduct(iProduct);
		
		const double vrj = reaction.GetStoichiometryOfProduct(iProduct);
		
		const double aj = GetActivityOfSpecies(idxPrimarySpecies[j]);
		
		Qr *= pow(aj, vrj);
	}
		
	return Qr;
}

const double System::KineticReactionQuotient(unsigned iKinetic) const // TODO: Modify the design of the class such that there is only one method to compute reaction quotient.
{
	const KineticReaction* reaction = kineticReactions[iKinetic];
	
	const unsigned r = reaction->GetIndexOfReactant();
	
	const double ar = GetActivityOfSpecies(idxKineticSpecies[r]);
	
	double Qr = 1.0/ar;
	
	for(unsigned iProduct = 0; iProduct < reaction->GetNumOfProducts(); ++iProduct)
	{
		const unsigned j = reaction->GetIndexOfProduct(iProduct);
		
		const double vrj = reaction->GetStoichiometryOfProduct(iProduct);
		
		const double aj = GetActivityOfSpecies(idxPrimarySpecies[j]);
		
		Qr *= pow(aj, vrj);
	}
		
	return Qr;
}



//
//const vector<double> System::GenerationRatesOfKineticSpecies(const vector<double>& molesOfMinerals) const
//{
//	vector<double> rates(GetNumOfKineticSpecies());
//	
//	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
//	{
//		const unsigned phaseIndexOfSpecies = GetPhaseIndexOfSpecies(GetIndexOfKineticSpecies(iKinetic));
//		
//		const double T = GetPhase(phaseIndexOfSpecies).GetTemperature();
//		
//		const double P = GetPhase(phaseIndexOfSpecies).GetPressure();
//		
//		rates[iKinetic] = kineticReactions[iKinetic].ReactionRate(T, P, molesOfMinerals[iKinetic], *this);
//	}
//	
//	return rates;
//}
//
//const vector<double> System::GenerationRatesOfPrimarySpecies(const vector<double>& ratesOfKineticSpecies) const
//{
//	vector<double> rates(GetNumOfPrimarySpecies(), 0.0);
//	
//	for(unsigned iKinetic = 0; iKinetic < GetNumOfKineticSpecies(); ++iKinetic)
//	{
//		const MineralReaction& mineralReaction = GetKineticReaction(iKinetic);
//		
//		for(unsigned iProduct = 0; iProduct < mineralReaction.GetNumOfProducts(); ++iProduct)
//		{
//			const unsigned iPrimary = mineralReaction.GetIndexOfProduct(iProduct);
//			
//			const double vkj = mineralReaction.GetStoichiometryOfProduct(iProduct);
//			
//			rates[iPrimary] -= vkj * ratesOfKineticSpecies[iKinetic];
//		}
//	}
//	
//	return rates;
//}
